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SUMMARY 


A study was conducted to develop an implicit method for integrating the equations of 
motion of a lumped-mass model of a rotor bearing system. The approach was, first, to use 
a Nordsieck-like numerical integration directly on the second-order equations of motion 
and, second, to assume that the forces and torques on the rotor are functions of the position 
and velocity at the point of application and its nearest axial neighbors. This allows the 
variables to be arranged so that the Jacobian of the set of nonlinear equations is block 
tridiagonal. Therefore the computational time is proportional to the number of elements in 
the rotor dynamics model rather than to the cube of the number. Numerical stability was 
demonstrated for any linearized homogeneous mode. 

To decrease computational time, a closed-form solution to the short-bearing theory 
was derived for a damper with arbitrary motion. Explicit results were presented for no 
cavitation and for full cavitation. 

The vast amount of data generated by the computer code was displayed in a motion 
picture showing an oblique view of the rotor bearing system. The motion of the rotor could 
be easily interpreted. 

An example problem of a rotor accelerating through three critical speeds with 19 
mass elements in the rotor dynamics model took 0.7 second of central processing unit time 
per time step on an IBM 360-67 computer in a time-sharing mode. The mode shapes at the 
first and third critical speeds were similar to the predicted mode shapes and occurred at 
the predicted speed. Because of the unbalance distribution, the second mode was not 
excited. Above the third critical speed the rotor bearing system operated as a self- 
centering device. This was also observed experimentally. 

The computer code, for the first time, allows us to look at a complex rotor bearing 
system with nonlinear transients and displays the vast amount of results in an easily 
understood motion-picture format. A 10-minute 16-millimeter, color, sound motion-picture 
supplement is available on loan. 


INTRODUCTION 

Nonlinear transients that are important in flexible, rotating equipment are difficult 
to analyze. Such things as blade tip rubs, spline friction, and squeeze-film dampers are 
difficult to predict with a linear model. Some of the transients that are important are 
locked rotor starts, blade loss, and rapid deceleration due to bearing failures. 

There are two basic methods for studying transient rotor oynamics. The first method 
is the modal method (refs. 1 and 2). It is best suited to linear rotor bearing systems running 
at a constant speed. The second method is the direct integration of the equations of mo- 
tion. It can be applied easily to nonlinear systems that are varying' in speed. The problem 
with the direct method is that it is limited by either computer running time or numerical 
stability. 

The equations of motion for rotor dynamics can be integrated directly in either of 
two ways, explicit or implicit integration. The explicit integration method solves the 
equations of motion at the present time for higher order derivatives and then extrapolates 
the displacements and velocities with a Taylor series to the advanced time (ref. 3). The 
implicit method solves the equations of motion (implicitly) at the advanced time step for 



the displacements and velocities, such that an extrapolation backward in time gives the 
previous results. 

The explicit method tends to be unstable when the product of the critical frequency 
(for any mode numerically possible! and the time step is large (ref. 4). Since the highest 
frequency is related to the square of the number of elements in the rotor dynamics moael, 
the computational time will be related to the square of the number of elements. Ap- 
proximately five or six elements seems to be a practical limit to the explicit method (ref. 2); 
that is, it can only be applied to simple assemblies. 

In contrast to the explicit method, the implicit method tends to be stable for large 
time steps (ref. 51; but it requires the solution of a large number of nonlinear simultaneous 
equations at each time step. For every element in the rotor dynamics model there are four 
degrees of freedom. For each degree of freedom there is an associated displacement and 
velocity. Therefore the total number of nonlinear equations to be solved at each time step 
is eight times the number of elements in the rotor dynamics model. The number of com- 
putations necessary to solve these equations is proportional to the cube of the number of 
equations. Therefore the computing time is proportional to the cube of the number of 
elements. 

This study was conducted to develop an implicit method for integrating the equations 
of motion in a reasonable amount of computing time. The approach is, first, to use a 
N ordsieck-like numerical integration directly on the second-order equations of motion* 
and, second, to assume that the forces and torques on the rotor are functions of the position 
and velocity of the point where the force or torque is applied and its nearest axial 
neighbors. This allows the variables to be arranged so that the Jacobian of the set of 
nonlinear equations is block tridiagonal. The computing time is proportional to the number 
of elements in the rotor dynamics model rather than to the cube of the number of elements. 

Besides the problems associated with integrating the equations of motion, there is a 
problem of describing the nonlinear damper force at each instant of time for an arbitrary 
orbit. In the past this was done by numerically integrating the Reynolds equation around the 
damper (ref. 6). This required a considerable amount of computing time. As an aside, a 
closed-form solution to the short-bearing theory was derived for a damper with arbitrary 
motion. 


SYMBOLS 

A coefficients used in partial-fraction expansion 

a property of shaft between mass stations defined in eq. (18a) 

b property of shaft between mass stations defined in eq. (18b) 

C radial clearance 

c property of shaft between mass stations defined in eq. (18c) 

D diameter 


iThis method of numerical integration was developed by Frank J. Zeleznik of the 
Lewis Research Center. For a set of first-order equations, it reduces to Gear's 
method. 
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modulus of elasticity 

force 

torque 

clearance in direction n 
moment of inertia 
index 
index 

axial length between mass stations 
mass of a rotor segment 
number of rotor segments 
radial direction at angle 6 
order of error in Taylor series 
pressure 

order of Taylor series 
radial displacement 
stability matrix 
element of S 
time 

time step 
defined in eq. (5) 

nondimensional velocity of journal in rotating coordinates 

real part of radial displacement 

imaginary part of radial displacement 

independent variable 

axial coordinate 

given set of constants 
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r angle defined in eq. (431 

e eccentricity ratio 

£ damping ratio 

6 circumferential angle 

A eigenvalue of stability matrix 

p viscosity 

co frequency 

Subscripts: 

B bearing 

J journal 

P polar 

T transverse 

+ associated with nearest axial neighbor or root of eq. (451 

0 start of integration 

1 end of integration 
Superscripts: 

(’1 time derivative 

( 1' axial derivative 

(") average or conjugate 

(k) kth time derivative 

-*■ vector 

unit vector 

NUMERICAL INTEGRATION 

Given an arbitrary function Z k (tl whose derivatives exist, Z^(t), a 1 'ay lor series 
expansion can be written: 


4 



Z k (t + At) = 


( 1 ) 


q-k 


(At) j fj) 


V ■ ■ Z' J ' (t) + 0 

j I k q-k 

J-0 

with Lagrange's remainder of order Oq-^. If the arbitrary function is chosen as 


7 = (At) K (k) 

^k k! r 


( 2 ) 


the Taylor series for this function becomes 


Z k (t + At) 


■ i (:>» 

j-0 


(t) + 0, 


(3) 


where the binomial coefficients are defined as 


^ _ J k! (j - k) : 


for j > k 
for j < k 


(4) 


If the form of the remainder is chosen as 


0 = a u 

q k 


the Taylor series becomes 


(5) 


Z k (t + At) 


■ 2 (!)>. 
J-o 


(t) + a k u 


( 6 ) 


where is a given set of constants and u can be determined from the equation of motion 
at the advanced time. The equation of motion at the advanced time is 

£F(r, r, r, t + At) = 0 ( 7 ) 

From the definition of z, the various derivatives become 


r 


(k) „ 


k! 

(At) k 


( 8 ) 
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Substituting for the various derivatives into the equation of motion and knowing the values 
at the previous time result in the equation of motion being a function of 

X)f(u, t + At) - 0 (9) 


This equation can be solved for u and, from this value of u, the remainder can be used as 
an error estimate to control the time step. 


N UMEitICAJL STABILITY 

The analysis of the stability of the numerical integration technique assumes a model 
of a rotor bearing system that is linearized at some instant of time. The homogeneous 
equation of motion for any mode is 


r + 2oo£r + a; r = 0 


( 10 ) 


where u> is the natural frequency and £ is the damping ratio for the mode. For every 
mode that is numerically possible, with nonnegative damping ratio, the amplitude must 
either remain constant or decay in time. The numerical integration is defined as unstable 
if the amplitude grows in time. 

From the definition of z the modal equation becomes 


2Z 2 + 2to At £ Z 1 + (u) At) Z = 0 (11) 

Substituting the Taylor series into the modal equation at the advanced time results in 

. 2 " 


-ze 


(j - 1) + 2,1(0 At g + (to At)' 


2c» 2 + 2a^w At £ + Oq((o At)' 


z j(t) 


( 12 ) 


Jj = 0 


For this value of u, the Taylor series expresses the solution at the advanced time in terms 
of the solution at the present time as 


Z k (t + At) 


-I ftp 


j=0 


j (j - 1) + 2j(0 At ? + (to At) 


2a 2 + 2a^io At £+ a Q ((o At)‘ 



(13) 


Defining the matrix element sjg to be 


kj 


(j\_^ 

j (j - 1) + 2j(0 At £+ (to At) 

w 

2a 2 + 2a^to At £+ oig((o At)^J 


(14) 
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and the q-dimension vector Z gives the eigenvalue equation as 

SZ = XZ (15) 

If the | X | > 1, the amplitude grows ana the method is numerically unstable. For q = 2, 
the given ot's are ag =2, oq = 3, and a2 = 1. In the limit as ui At ■+ the maximum 
| X | ■+ 0. Therefore if the time step At is much larger than w“l for a mode, the 
amplitude of that mode will approach zero. If a mode is to have a nonzero amplitude, uh 
must be small. In the limit as to At 0, the maximum |X| -*■ 1. Therefore the method 
is numerically stable in the two limits. 


EQUATIONS OF MOTION 

A model of the shaft showing the complex number representation of the radial 
displacement r is shown in figure 1. The radial displacement is the distance between the 
shaft centerline and the axis of rotation. It can be represented by 


r = x + iy (16a) 

where the real and imaginary axes are fixed in space perpendicular to the axis of rotation. 
The slope of the shaft along the axis of rotation is 

r' = x' + iy' (16b) 

The position of the shaft is then described by r and r' at all the axial locations. 

The lumped-mass model of a rotor divides the rotor into N segments. The mass and 
inertia of each segment are assumed to be concentrated at a point. These points are then 
assumed to be connected by massless elastic beams that model the stiffness of the rotor. 

The equations of motion for the lumped-mass model were derived in reference 7. The 
sum of the forces y,F at a point, must be zero, where 

T.F = -mr + a_r_ - (a_ + a + )r + a + r + 

+ b r' + (b - b . )r ' - b r! + F = 0 (17a) 

“ - “ + + + 

and the sum of the torques y^G about a point must be zero, where 


= -I T r’ + iul p r ' - uil p r’ 

- b_r_ + (b_ - b + )r + b + r + 

- c_v'_ - 2(c_ + c + )r ' - c + r| + G (17b) 
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The + or - refer to the next or previous axial location; and a, b, and c are properties of 
the shaft between these locations: 


a = 12EI/L 3 

(18a) 

b = 6EI/L 2 

(18b) 

c = 2EI/L 

(18c) 


If the nonlinear force F and the nonlinear torque G are functions of displacements 
and velocities of the point and its nearest neighbors, the y^F and Y'p are functions of the 
displacements and velocities of the point and its nearest neighbors. If the Taylor series of 
the numerical integration technique is substituted into the equations of motion for the 
acceleration, velocity, and displacement, the form of the equations of motion becomes 


=2>(u_. 

u > u +» u l> u '» u +) 

(19a) 

ii 

M 

O 

/•— N 

c 

u > u +> u l» u '» u +) 

(19b) 


These equations form a set of 2N complex nonlinear equations in 2N unknowns. 
These equations are solved by rewriting them as 4N real equations in 4N real unknowns 
and then using a Newton-ftaphson iterative technique to obtain a numerical solution. The 
N ewton-Raphson technique assumes a solution, linearizes the equations about that solution, 
and then solves the linear set of equations for a correction to the assumed solution. The 
form of the equations of motion results in the linear set of equations being block tri- 
diagonal (fig. 2). The block-tridi agonal form allows the set of equations to be solved in a 
very efficient manner. The computing time is proportional to N rather than to as in 
the general method. 


SQUEEZE-FILM DAMPER BEARING 

The configuration of the squeeze-film damper is shown in figure 3. The same con- 
figuration can be used to analyze journal bearings where the journal and the bearing are 
allowed to rotate. If coj is the rotational speed of the journal and tog is the rotational 
speed of the bearing, the average rotational speed is 

0) _ _1 (ojg + Wj) (20) 

For a damper this average rotational speed would be zero. 

If C is the radial clearance and T is the displacement of the journal center with 
respect to the bearing center, the clearance h in the direction n is 

■f A 

h = c - r • n (21) 
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( 22 ) 


If r is the velocity of the journal, 


If n is at an angle 0, 


3h 4- '• 

- = - r S n 


3n 

30 = k X n 


(23) 


where k is a unit vector along the axis of the damper bearing in the z-direction. If co is 
defined as 


u = uk 


(24) 


then 


0 ) -^ = (w X r) • n 


(25) 


The Reynolds equation for the short, plain damper journal bearing is presented in 
reference 6 as 


3_ / h 3 3 p\ = - 3h 3h 

3z \12y 9Z J “ 36 3t 


(26) 


If the boundary conditions in the damper bearing are 

HO, 0, t) = 0 


(27a) 


P(L, 0, t) = 0 


(27b) 


and if h is not a function of z, 


P = - 


6pz(L - 2 ) .-*■ T 


(u> x r - r) • n 


(28) 


The eccentricity ratio is 


so that 


-*■ r 
E = C 


_ r 
E C 


(29) 

(30) 
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If V is defined as 


$ = 


*■>* 

a) x e 


( 31 ) 


P 


6 V 2 ( L ~ z > $ . n 
C 2 (l - e • n) 


(32) 


a — r 

The pressure is zero when n is perpendicular to V, and the pressure is greater than zero 
when V • n > 0. 

The force on the journal due to the pressure in a segment of the film extending from 
00 to 6 } is 


sr 

e 0 o 


D re-, rh 

F = - — / 1/ Pn dz d0 


(33) 


This expression for the force can be integrated axially and becomes 


_ uDL 3 S'® 1 (V • n)n 

2 J /- a . 3 

2C (1 - e • n) 

9 0 


de 


(34) 


The angular integral can be integrated by transforming the integral to the complex 
plane (fig. 4). Let V be in the real direction, e be at an angle <f>, and n be at an angle 
6 so that 


V = I V [ 

(35a) 

e = [ e | e i( ^ 

(35b) 

id 

n = e 

Differentiating the expression for n yields 

(35c) 

d0 = -in ^ dn 

(36) 

and using the definition of the complex cosine yields 


^ - V(n + n 1 ) 

V • n = — 1 

(37) 


2 


n = 


(ne + en ^) 


(38) 
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The expression for the force becomes 


n 2 (n 2 + l)dn 


F = 



2n + e) 3 


(39) 


where the integral is around the unit circle from nQ to nj, where 


n 


0 



n 


1 



(40a) 

(40b) 


The pressure is zero at n = + i, and the pressure is greater than zero when /<?g(n) is greater 
than zero. 

For no cavitation the integral extends completely around the journal; and by using the 
theory of residues, the force becomes 


F = 



(41) 


For cavitation the integral extends from -i to +i; and by using a partial-fraction tech- 
nique, the force becomes 


F 



(1 + njV 1 


(-n ) j 2 A 


+ 


JU. 


(1 + n 2 ) j 1 


J 


(42) 


where I is defined as 


r 


tan 


-1 


Vi - ui 2 

' -&(e) 


(43) 


and r is in the first or second quadrant. The partial-fraction expansion coefficients are 


A 


±3 


n 2 (n 2 + 1) 
(n± - nJ> 3 


(44a) 
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n 


n_ 

+ 


(44b) 


2n ± (2n^ + 1) 

S±2 = ~ 73 ~ 

(n ± - n_) 


3A 


±3 


6n^ + 1 

Sfcl = 7 

(n± - n_) 


3A 


+2 


3A 


±3 


n , - n_ , .2 

± + (n + - n_) 


(44c) 


where the roots of the denominator of the force equation are 


n+ = 



(45) 


DISCUSSION OF EXAMPLE 

The rotor bearing- system described in reference 8 was used as the example problem. 
This rotor bearing system consisted of a shaft with three disks mounted on two axially pre- 
loaded ball bearings (fig. 5). The bearings were mounted in a squeeze-film damper journal, 
and the journal had a centering spring. 

The first three critical speeds for the rotor bearing system without oil in the damper 
are shown in figure 6. All the modes are bent-shaft modes. The "classical" hierarchy only 
applies to stiff shafts; therefore, the classical mode shapes do not characterize the actual 
mode shapes. The first mode, about 7581 rpm, classically would be the cylindrical mode. 
But in this case, it has a large amount of bending outward near the shaft center. The sec- 
ond mode, about 9235 rpm, classically would be the conical mode. In this case, it has a 
slight amount of bending outward near the shaft ends. The third mode, about 11 248 rpm, 
classically would be the bending mode. In this case, it has a large amount of bending 
throughout the shaft. 

Experimentally the rotor was accelerated from rest through the three critical 
speeds. The Lissajous patterns for the three disks were displayed on three side-by-sioe 
cathode ray tubes. A motion picture was taken of the CRT's plus a speed counter. The 
Lissajous patterns at the three critical speeds are shown on figures 7 to 9. The three 
critical speeds occurred at about the predicted speeds, and the Lissajous patterns corre- 
sponded to the three mode shapes. 

The rotor bearing system was modeled by using 19 elements. The rotor was assumed 
to have a uniform, in-line unbalance, with a mass eccentricity of 0.00254 centimeter (1 
mil). The equations of motion for this system were programed in FORTRAN IV on an IBM 
360-67 computer in a time-sharing mode. The equations of motion were directly integrated 
by the implicit integration method, with a fixed time step of 0.12 millisecond. The tran- 
sient analyzed was the rotor accelerated from rest through the three critical speeds. The 
rate of acceleration was 8727 rad/sec^. Each time step took about 0.7 second of CPU time. 

The output at each time step of the calculation was displayed on a CRT. The display 
showed an oblique view of the rotor bearing system, with the bearing centerline as the ob- 
lique axis. The transverse vibration is indicated by a series of dots. Each dot represents a 
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location of an element in the rotor dynamics model. The scale of the transverse vibration 
exaggerates the amplitude of the vibration. The display on the CRT was photographed at 
each time step. These photographs were then shown as a motion picture. 

The computer-generated displays on the CRT at the first and third critical speeds and 
at a speed much greater than the third critical speed are shown in figures 10 to 12. The 
mode shapes at the first and third critical speeds were similar to the predicted mode shapes 
and occurred at the predicted speed. Because of the unbalance distribution, the second 
mode was not excited. The only indication of the second mode was a traveling wave super- 
imposed on the first mode shape. This traveling wave decayed when the rotor went through 
the third critical speed. Above the third critical speed, the rotor bearing system operated 
as a self-centering device. The mass centerline coincided with the bearing centerline. 
Therefore, the rotor displacement was uniform and in line, with an amplitude of 0.00254 
centimeter (1 mil). 

In conclusion, this computer code for the first time allows us to look at complex rotor 
bearing systems experiencing nonlinear transients and displays the vast amount of results in 
an easily understood motion-picture format. A 10-minute, 16-millimeter, color, sound 
motion-picture supplement is available, on loan, that shows the test data and the computer- 
made motion picture. 


CONCLUSIONS 

An implicit method for integrating the equations of motion for a lumped-mass model 
of a rotor dynamic system was developed. The following conclusions were drawn: 

1. The method was numerically stable for any time step. 

2. An error estimate was available to control the size of the time step. 

3. The computational time was proportional to the number of elements in the rotor 
dynamics model rather than to the cube of the number. 

4. An example problem with 19 mass elements in the rotor dynamics model took 0.7 
second of central processing unit time per time step an an IBM 360-6 7 computer in a 
time-sharing mode. 

For the first time, this code allows the simulation of a complex rotor bearing system 
experiencing nonlinear transients and displays the vast amount of results in an easily 
understood motion-picture format. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 19, 1979, 

505-04. 
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Figure 1. - Model of shaft showing complex number representation 
of radial displacements. 



Figure 2. - Newton-Raphson technique that leads to a linear set of block-tridiagonal equations. 
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Data acquisition 



Figure 5. - Rotor bearing system used as example problem. 
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Figure 6. - Undamped critical speeds. 













I 



Figure 9. - Rotor passing through third critical speed. 
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Figure 10. - Oblique view of rotor centerline as rotor 
passes through first critical speed. 
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Figure 11. - Oblique view of rotor centerline as rotor 
passes through third critical speed. 
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Figure 12. - Oblique view of rotor centerline at rotor 
speeds much greater than third critical. 
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